{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "<h1><font color=\"darkblue\">Robust Statistics</font></h1>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Datasets\n",
    "\n",
    "- For example, a set of $n$ scalar measurements \n",
    "\n",
    "> $\\displaystyle D = \\big\\{x_i: \\ i=1\\dots{}n\\big\\} $"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How to characterize the data?\n",
    "- Location\n",
    "- Dispersion\n",
    "- Shape?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Maximum Likelihood Estimation\n",
    "\n",
    "- MLE looks for the optimal parameter of the likelihood function\n",
    "\n",
    "> Worked well before. Using Gaussian likelihood, the MLE for the location was (weighted) average of measurements. We also used the sample variance for dispersion.\n",
    "\n",
    "- Problems with outliers when using, e.g., Gaussian likelihood\n",
    "\n",
    "> Just a few outliers can throw these estimates big time!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Location\n",
    "\n",
    "- MLE looks for the parameter for the optimial likelihood function\n",
    "\n",
    ">$\\displaystyle L(\\mu) = p(D \\lvert \\mu) = \\prod_{i=1}^n \\ell_i(\\mu)$\n",
    "\n",
    "> with $\\quad\\ell_i(\\mu) = f_0(x_i-\\mu)$\n",
    "\n",
    ">$\\displaystyle \\hat{\\mu} = \\arg \\max_{\\mu} \\prod_i\\,f_0(x_i-\\mu)$\n",
    "\n",
    "- Very general but for a special family of functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### How to solve it?\n",
    "\n",
    "- Solve it as before: take its log and differentiate\n",
    "\n",
    ">$\\displaystyle \\hat{\\mu} = \\arg \\min_{\\mu} \\sum_{i=1}^n \\rho(x_i-\\mu)$\n",
    "\n",
    "> where $\\quad\\rho = -\\log f_0$\n",
    "\n",
    "- At the minimum the derivative vanishes\n",
    "\n",
    ">$\\displaystyle \\sum_{i=1}^n \\rho'(x_i-\\hat{\\mu})=0$\n",
    "\n",
    "- And that's it! Almost..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What $\\rho$?\n",
    "\n",
    "- For Gaussian likelihood\n",
    "\n",
    ">$\\displaystyle \\rho(t) = \\frac{t^2}{2}$\n",
    "\n",
    "> and we have\n",
    "\n",
    ">$\\displaystyle \\sum_{i=1}^n (x_i-\\hat{\\mu})=0$\n",
    "\n",
    "- So $\\hat{\\mu}$ is the **mean** (as before)\n",
    "\n",
    "> What if $\\rho(t)=t^2$ instead?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Robust $\\rho$-function\n",
    "\n",
    "- Less sensitive to outliers, e.g.,\n",
    "\n",
    ">$\\rho(t) = \\lvert\\,t\\,\\rvert$\n",
    "\n",
    "> and we have\n",
    "\n",
    ">$\\displaystyle \\sum_{i=1}^n \\textrm{sgn}(x_i-\\hat{\\mu})=0$\n",
    "\n",
    "- So $\\hat{\\mu}$ is the **median** \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Different $\\rho$-functions\n",
    "\n",
    "- Many possible functions to choose from, e.g., \n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/c/c1/RhoFunctions.png\" width=400 align=left >\n",
    "\n",
    "\n",
    "<!-- <img src=\"files/RhoFunctions.png\" width=400 align=left \\> -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The trick\n",
    "\n",
    "- Define a new function $W$ \n",
    "\n",
    ">$\\displaystyle W(t) = \\frac{\\rho'(t)}{t}$\n",
    "\n",
    "> so we have\n",
    "\n",
    ">$\\displaystyle \\sum_{i=1}^n W(x_i\\!-\\!\\hat{\\mu})\\,(x_i\\!-\\!\\hat{\\mu})=0$\n",
    "\n",
    "- If we had constant $w_i$ weights\n",
    "\n",
    ">$\\displaystyle \\sum_{i=1}^n w_i\\,(x_i\\!-\\!\\hat{\\mu})=0$\n",
    "\n",
    "> the solution would be the weighted average\n",
    "\n",
    ">$\\displaystyle \\hat{\\mu} = \\frac{\\sum w_ix_i}{\\sum w_i}$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Iterative method\n",
    "\n",
    "- Intuitive and efficient\n",
    "\n",
    ">0. Obtain initial estimate of $\\hat{\\mu}$, e.g., median\n",
    "0. Assign $w_i=W(x_i\\!-\\!\\hat{\\mu})$ weights to the measurements\n",
    "0. Calculate the weighted average\n",
    "0. Repeat steps 2 and 3 until convergence\n",
    "\n",
    "- Very fast convergence in practice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Wc(t,c=1): # Cauchy\n",
    "    return 1 / (1 + np.square(t/c))\n",
    "\n",
    "def Wh(t,k=1): # Huber\n",
    "    w = np.ones_like(t)\n",
    "    abst = np.abs(t)\n",
    "    i = (abst > k)\n",
    "    w[i] = k / abst[i]\n",
    "    return w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "t = np.linspace(-3,3,200); plt.ylim(0, 1.1);\n",
    "[plt.plot(t,w) for w in [Wc(t), Wh(t), Wh(t,0.5)]];"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.8911863410757943, -0.10991278679212746)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.random.seed(42)\n",
    "x = np.random.randn(100,1)\n",
    "x[0] = 200 # outlier\n",
    "x.mean(), x[1:].mean()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2.54789391e-05 1.95363252e-01 3.92728740e-01]]\n",
      "[[2.51313636e-05 6.93908140e-01 9.85384110e-01]]\n",
      "[[2.50406264e-05 9.15873333e-01 8.10917106e-01]]\n",
      "[[2.50028007e-05 9.77427536e-01 7.13297811e-01]]\n",
      "[[2.49867618e-05 9.92351907e-01 6.72299766e-01]]\n",
      "[[2.49799665e-05 9.96343371e-01 6.55237403e-01]]\n",
      "[[2.49770931e-05 9.97597656e-01 6.48091352e-01]]\n",
      "[[2.49758796e-05 9.98049158e-01 6.45086280e-01]]\n",
      "[[2.49753673e-05 9.98225747e-01 6.43820136e-01]]\n",
      "[[2.49751511e-05 9.98297776e-01 6.43286210e-01]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-0.09697112])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAGiFJREFUeJzt3XtwVGWaBvDnTYckBOUWAmJCEmoGEbkJRASVkOFy0HHVUhxFIxbqTobFkZn1Ppsad2ucuDslirBaQBRk0YxowZSjonIzEFRAwkUJoIJKQrhICHKRkGCSd//oJgbSSXdId74+5zy/qq50n/66+7ElT598p/trUVUQEZGzRJkOQEREocdyJyJyIJY7EZEDsdyJiByI5U5E5EAsdyIiB2K5ExE5EMudiMiBWO5ERA4UbeqBu3XrpmlpaaYenojIljZv3nxEVRMDjTNW7mlpaSgqKjL18EREtiQiJcGM47QMEZEDsdyJiByI5U5E5EAsdyIiB2K5ExE5EMudiMiBWO5ERA7EciciciCWOxGRA7HciYgciOVORORALHciIgdiuRMROVDAcheRXiJSICI7RWSHiPzBzxgRkdkiskdEvhCRoeGJS0REwQhmyd8aAI+o6hYRuRjAZhFZqao7G4y5AUAf3+lqAHN8P4mIyICAe+6qelBVt/jOnwSwC0DSecNuAbBIvTYA6CwiPUOeloiIgtKiOXcRSQMwBMDG865KArCvweUyNH4BgIhki0iRiBSVl5e3LCkREQUt6HIXkYsALAXwR1U9cSEPpqp5qpququmJiQG/JYqIiC5QUOUuIu3gLfZ8Vf2HnyH7AfRqcDnZt42IiAwI5t0yAmA+gF2q+nwTw94BcK/vXTMjABxX1YMhzElERC0QzLtlrgUwGcB2Ednm2/YfAFIAQFXnAngfwK8B7AFQCeC+0EclIqJgBSx3Vf0YgAQYowAeDFUoIiJqHX5ClYjIgWxV7vn5+UhLS0NUVBTS0tKQn59vOhIRUUQKZs49IuTn5yM7OxuVlZUAgJKSEmRnZwMAsrKyTEYjIoo4ttlzz8nJqS/2syorK5GTk2MoERFR5LJNuZeWlrZoOxGRm9mm3FNSUlq0nYjIzWxT7rm5uYiPjz9nW3x8PHJzcw0lIiKKXLYp96ysLOTl5SE1NbV+28yZM3kwlYjID9uUO+At+L179+LTTz8FAHTp0sVwIiKiyGSrcj/rqquuQqdOnbBixQrTUYiIIpItyz06Ohpjx47FihUr4F35gIiIGrJluQOAZVkoLS3F119/bToKEVHEsXW5A+DUDBGRH7Yt9969e+OXv/wly52IyA/bljvg3XsvKCjAmTNnTEchIoooti/3U6dOYf369aajEBFFFFuXe2ZmJjweD6dmiIjOY+ty79SpE0aMGMFyJyI6j63LHfBOzWzevBlHjhwxHYWIKGI4otxVFatXrzYdhYgoYti+3NPT09G5c2dOzRARNWD7cudSBEREjdm+3AHv1ExZWRm+/PJL01GIiCKCI8p9/PjxALgUARHRWY4o9969e6NPnz4sdyIiH0eUO+CdmlmzZg2qq6tNRyEiMs5R5V5ZWVn/LU1ERG7mmHLPzMxEdHQ0p2aIiOCgcu/YsSNGjhzJcicigoPKHfBOzWzZsgXl5eWmoxARGeW4cgeAVatWGU5CRGSWo8p92LBh6NKlC6dmiMj1HFXuHo8H48aN41IEROR6jip3wDs1c+DAAezcudN0FCIiYxxX7lyKgIjIgeWempqKvn37styJyNUClruILBCRwyJS3MT1mSJyXES2+U5PhT5my1iWhbVr16Kqqsp0FCIiI4LZc18I4PoAY9ap6pW+019aH6t1LMvC6dOn8cknn5iOQkRkRMByV9VCAEfbIEvIZGZmol27dpyaISLXCtWc+0gR+VxEPhCR/iG6zwt20UUX4ZprrmG5E5FrhaLctwBIVdXBAP4XwNtNDRSRbBEpEpGicC8RYFkWtm3bhu+//z6sj0NEFIlaXe6qekJVf/Sdfx9AOxHp1sTYPFVNV9X0xMTE1j50s7gUARG5WavLXUQuERHxnR/uu8+K1t5vaw0ZMgQJCQmcmiEiV4oONEBE3gCQCaCbiJQB+E8A7QBAVecCuB3Av4lIDYDTACZpBHz2//ylCHyvP0RErhCw3FX1rgDXvwjgxZAlCiHLsvDmm2+iuLgYAwcONB2HiKjNOO4Tqg1xKQIicitHl3uvXr3Qr18/ljsRuY6jyx3wTs0UFhbi9OnTpqMQEbUZV5R7VVUVPv74Y9NRiIjajOPLffTo0VyKgIhcx/Hl3qFDB1x33XUsdyJyFceXO+Cdmvniiy9w8OBB01GIiNqEa8od4FIEROQerij3K6+8Et26dePUDBG5hivKPSoqCuPHj8fKlStRV1dnOg4RUdi5otwB79TM999/j+3bt5uOQkQUdq4pdy5FQERu4ppyT0pKQv/+/VnuROQKril3wDs1s27dOlRWVpqOQkQUVq4r9+rqaqxbt850FCKisHJVuWdkZCAmJoZTM0TkeK4q9/j4eIwaNYrlTkSO56pyB7xTM8XFxThw4IDpKEREYePKcgeAlStXGk5CRBQ+riv3QYMGoXv37pyaISJHc125cykCInID15U74J2aKS8vx+eff246ChFRWLiy3LkUARE5nSvLvWfPnhg4cCDLnYgcy5XlDninZj7++GOcOnXKdBQiopBzdbmfOXMGhYWFpqMQEYWca8t91KhRiI2N5dQMETmSa8u9ffv2yMjIYLkTkSO5ttwB79TMzp07UVZWZjoKEVFIub7cAWDVqlWGkxARhZary33gwIHo0aMHp2aIyHFcXe4iAsuyuBQBETmOq8sd8E7NHDlyBNu2bTMdhYgoZFxf7uPGjQPApQiIyFlcX+6XXHIJBg8ezHInIkdxfbkDXIqAiJwnYLmLyAIROSwixU1cLyIyW0T2iMgXIjI09DHDy7Is/PTTT1i7dq3pKEREIRHMnvtCANc3c/0NAPr4TtkA5rQ+Vtu67rrrEBcXx6kZInKMgOWuqoUAjjYz5BYAi9RrA4DOItIzVAHbQlxcHEaPHs1yJyLHCMWcexKAfQ0ul/m22YplWdi1axf27dsXeDARUYRr0wOqIpItIkUiUlReXt6WDx3Q2aUIVq5caTgJEVHrhaLc9wPo1eBysm9bI6qap6rpqpqemJgYgocOnf79+6Nnz56cmiEiRwhFub8D4F7fu2ZGADiuqgdDcL9tquFSBLW1tabjEBG1SjBvhXwDwHoAfUWkTEQeEJGpIjLVN+R9AN8C2APgZQDTwpY2zCzLwtGjR7F161bTUYiIWiU60ABVvSvA9QrgwZAlMqjhUgTp6emG0xARXTh+QrWB7t27Y8iQIZx3JyLbY7mfx7IsfPrppzh58qTpKEREF4zlfh4uRUBETsByP8+1116L9u3bc2qGiGyN5X6e2NhYZGZmstyJyNZY7n5YloWvvvoKJSUlpqMQEV0QlrsfXIqAiOyO5e5Hv379kJSUxKkZIrItlrsfZ5ciWLVqFZciICJbYrk3wbIs/PDDD9i8ebPpKERELcZyb8K4ceMgIpyaISJbYrk3oVu3bhg6dCjLnYhsieXeDMuysH79epw4ccJ0FCKiFmG5N8OyLNTU1GDNmjWmoxARtQjLvRkjR45Ehw4dODVDRLbDcm8GlyIgIrtiuQdgWRZ2796N7777znQUIqKgsdwD4FIERGRHLPcA+vbti169enFqhohsheUewNmlCFavXo2amhrTcYiIgsJyD4JlWTh27BiKiopMRyEiCgrLPQhjx47lUgREZCss9yAkJCQgPT2d5U5EtsFyD5JlWdiwYQOOHz9uOgoRUUAs9yBZloXa2loUFBSYjkJEFBDLPUgjRozgUgREZBss9yDFxMTgV7/6FcudiGyB5d4ClmXhm2++wTfffGM6ChFRs1juLcClCIjILljuLXDZZZchJSWFUzNEFPFY7i3ApQiIyC5Y7i1kWRZOnDiBzz77zHQUIqImsdxbiEsREJEdsNxbqGvXrrjqqqtY7kQU0VjuF8CyLGzcuBHHjh0zHYWIyC+W+wUQEdTV1aFr165IS0tDfn6+6UhEROdgubdQfn4+nnvuOQCAqqKkpATZ2dkseCKKKEGVu4hcLyJficgeEXnSz/VTRKRcRLb5Tv8a+qiRIScnB5WVledsq6ysRE5OjqFERESNRQcaICIeAC8BGA+gDMAmEXlHVXeeN/RNVf19GDJGlNLS0hZtJyIyIZg99+EA9qjqt6p6BsBiALeEN1bkSklJ8bs9OTm5jZMQETUtmHJPArCvweUy37bzTRSRL0RkiYj08ndHIpItIkUiUlReXn4Bcc3Lzc1FfHx8o+0ejwcVFRUGEhERNRaqA6rvAkhT1UEAVgL4P3+DVDVPVdNVNT0xMTFED922srKykJeXh9TUVIgIUlNT8cgjj+DAgQMYM2YMDh8+bDoiEVFQ5b4fQMM98WTftnqqWqGq1b6LrwAYFpp4kSkrKwt79+5FXV0d9u7dixkzZuC9997D7t27kZmZiQMHDpiOSEQuF0y5bwLQR0R6i0gMgEkA3mk4QER6Nrh4M4BdoYtoD+PHj8cHH3yA0tJSjB49Gvv27Qt8IyKiMAlY7qpaA+D3AJbDW9pvqeoOEfmLiNzsGzZdRHaIyOcApgOYEq7AkWz06NFYsWIFDh8+jIyMDHz33XemIxGRS4mqGnng9PR0LSoqMvLY4VZUVATLshAfH4+PPvoIl112melIROQQIrJZVdMDjeMnVMMgPT0dBQUFOHPmDDIyMrBjxw7TkYjIZVjuYTJ48GCsWbMGUVFRyMzMxLZt20xHIiIXYbmH0RVXXIG1a9ciLi4OY8aMwaZNm0xHIiKXYLmHWZ8+fVBYWIjOnTtj3Lhx+OSTT0xHIiIXYLm3gd69e6OwsBA9evTAhAkTsGbNGtORiMjhWO5tJDk5GWvXrkVqaipuuOEGLF++3HQkInIwlnsb6tmzJ9asWYO+ffvi5ptvxrvvvms6EhE5FMu9jSUmJuKjjz7CoEGDcNttt2HJkiWmIxGRA7HcDejatStWrVqF4cOH48477+S3OBFRyLHcDenUqROWL1+OjIwMTJ48GQsWLDAdiYgchOVu0EUXXYRly5Zh/PjxeOCBBzBnzhzTkYjIIVjuhsXHx+Of//wnbrrpJkybNg0zZ840HYmIHIDlHgHi4uKwZMkSTJw4EQ8//DCeeeYZ05GIyOZY7hEiJiYGixcvxt13342cnBw89dRTMLViJxHZX7TpAPSz6OhoLFq0CHFxcXj66adRVVWFv/3tbxAR09GIyGZY7hHG4/Hg5ZdfRmxsLJ599llUVVXhhRdeQFQU/8giouCx3CNQVFQUXnrpJcTFxWHmzJmoqqrC3LlzWfBEFDSWe4QSETz33HNo3749nnnmGVRXV2P+/PmIjub/MiIKjE0RwUQEubm5iIuLw1NPPYXq6mq89tpraNeuneloRBThWO428Oc//xlxcXF4/PHHUV1djcWLFyM2NtZ0LCKKYJzEtYnHHnsMs2fPxttvv41bb70Vp0+fNh2JiCIYy91GHnroIcybNw8ffvghbrrpJpw6dcp0JCKKUCx3m8nOzsbChQtRUFCAYcOGISUlBVFRUUhLS+PqkkRUj3PuNnTvvfdi06ZNePHFF+u3lZSUIDs7GwCQlZVlKhoRRQjuuduUv29xqqysRE5OjoE0RBRpWO42VVpa6nd7SUkJVq9ejbq6ujZORESRhOVuUykpKX63R0VFYdy4cbj88ssxY8YMHDlypI2TEVEkYLnbVG5uLuLj48/ZFh8fj/nz5+O1115Djx498NhjjyEpKQlZWVlYt24dV5kkchGWu01lZWUhLy8PqampEBGkpqYiLy8PU6ZMwT333IN169Zh+/btyM7OxrJly5CRkYEBAwZg9uzZOHbsmOn4RBRmYmpvLj09XYuKiow8ttucOnUKb775JubNm4fPPvsM7du3x6RJk/C73/0Ow4cP55LCRDYiIptVNT3QOO65u0CHDh1w//33Y+PGjdi8eTMmT56Mt956CyNGjMDQoUMxb948nDx50nRMIgohlrvLnC3zAwcOYM6cOVBVTJ06FZdeeimmTp2KrVu3mo5IRCHAcnepjh071pf5hg0bcPvtt2PRokUYOnQorr76arz66quorKw0HZOILhDL3eVEpL7M9+/fj1mzZuHkyZO4//77cemll2L69OnYsWOH6ZhE1EIsd6rXpUuX+jIvLCzEjTfeiHnz5mHAgAEYNWoU8vPzUVVVZTomEQWB5U6NiEh9me/fvx/PPvssDh06hHvuuQfJycl49NFHsXv37vrx+fn5SEtL4wJmRBEkqLdCisj1AGYB8AB4RVX/57zrYwEsAjAMQAWAO1V1b3P3ybdC2ktdXR0KCgowd+5cvP3226ipqcHYsWPRv39/vPLKK+fMz8fHxyMvL48LmBGFQbBvhQxY7iLiAfA1gPEAygBsAnCXqu5sMGYagEGqOlVEJgG4VVXvbO5+We72dejQISxYsAB5eXkoKSnxOyY1NRV79+5t22BELhDKch8J4L9UdYLv8p8AQFX/u8GY5b4x60UkGsAhAInazJ2z3O2vtrYW7dq1a3JZgyeffBIDBgzAgAEDcPnll/OrAYlCINhyD2Y99yQA+xpcLgNwdVNjVLVGRI4DSADAVasczOPxICUlxe/ee7t27TBjxgzU1NTUj+3Tp0992Z89/eIXv0B0NL9WgCjU2vS3SkSyAWQDTa9qSPaSm5uL7Oxsv3Puv/nNb7B7924UFxfXn7Zt24alS5fW7+3HxsaiX79+jUo/JSWFyyIQtQKnZajV8vPzkZOTg9LSUqSkpCA3N7fZg6mVlZXYtWvXOaVfXFyMsrKy+jEXX3wx+vfv36j0u3fv3mTptzQHkR2Fcs49Gt4DqmMB7If3gOrdqrqjwZgHAQxscED1NlW9o7n7ZbnT+Y4dO4YdO3acU/jbt29HRUVF/Zhu3bo1Kvz+/ftj2bJlTf4FwYInJwlZufvu7NcAXoD3rZALVDVXRP4CoEhV3xGROACvARgC4CiASar6bXP3yXKnYKgqDh8+3Ggvv7i4GD/++GP9OI/Hg9ra2ka379GjB9auXYuEhAR06dIFHo+nLeMThVxIyz0cWO7UGqqK0tLS+qJ/8skng7pd586dkZCQgISEBHTt2rXReX/bOnbs2KL5f04PUTix3MlV0tLS/L5rp3v37nj++edx9OhRVFRUoKKiwu/548ePN3nfHo/Hb+n7O79x40Y8/fTTOH36dP3tTUwPRcoLDHOEPgfLnVwlPz+/VXPuNTU1OHr0aMAXgfPPNyzx5kRFRSE5ORmxsbGIi4sLeAp2nL+x7733Hh5++OGIeIGJhOMgTsvBcifXMbF3dvr06fqyP3r0KMaMGdPkh7qmTJmCqqoqv6fq6mq/20P5+yki6NixIzweD6Kjo+HxeII635KxDc///e9/P+e4yFkdO3bEtGnTEBUVBRGBiNSfD3ZbS27zxBNPnHNQ/qyEhATMmjWrfsqt4c9QbWt4/r777kN5eXmjHC39NDfLnciApqaHLmQ5BlXFTz/9FPAF4PwXid/+9rdN3uf06dNRW1uL2tpa1NTUBDwf7Dh/tzl06FCTOWJiYqCqqKurq//pViLSov/+UH5CNaJ8+CHQzL8ZIqMmTHgDCxcuxJkz1fXbYmJiMWHCFCxc2NJ7EwAxvlPHoG4RHQ0kJHyJiorGHw5PSOiGIUNmtDTEBXv00UebzDFjhv8cqtrkyXt9HVQbjwMUdXXeMWeL8ux1f/3rX3Hs2A+NHqtz5y544okn6sc2zNDUNuDsdWj0V1XDnP5u88ILs3DixNljO4cALAcQvg902q7ciSLZyJEjAQBLly5FRUUFEhISMHHixPrtbWHixIl+X2AmTpzYZhkuNEfDqYxQueOOO/zmuOOOO3DJJZeE9LGaM2nSpEY54uPjkZubG54HbO6VMpynYcOGKRGFx+uvv66pqakqIpqamqqvv/46czgkB7yfLwrYsZxzJyKykWDn3PlNTEREDsRyJyJyIJY7EZEDsdyJiByI5U5E5EDG3i0jIuUA/H+7cmDdwK/wa4jPx7n4fPyMz8W5nPB8pKpqYqBBxsq9NUSkKJi3ArkFn49z8fn4GZ+Lc7np+eC0DBGRA7HciYgcyK7lnmc6QITh83EuPh8/43NxLtc8H7accycioubZdc+diIiaYbtyF5HrReQrEdkjIsF9K7JDiUgvESkQkZ0iskNE/mA6k2ki4hGRrSLynukspolIZxFZIiJfisguEWm7dYcjjIj8u+93pFhE3hCRONOZws1W5S4iHgAvAbgBwBUA7hKRK8ymMqoGwCOqegWAEQAedPnzAQB/ALDLdIgIMQvAh6p6OYDBcOnzIiJJAKYDSFfVAQA8ACaZTRV+tip3AMMB7FHVb1X1DIDFAG4xnMkYVT2oqlt850/C+8ubZDaVOSKSDOBGAK+YzmKaiHQCkAFgPgCo6hlVPWY2lVHRANqLSDSAeAAHDOcJO7uVexKAfQ0ul8HFZdaQiKQBGAJgo9kkRr0A4HEA7v1Czp/1BlAO4FXfNNUrItLBdCgTVHU/gBkASgEcBHBcVVeYTRV+dit38kNELgKwFMAfVfWE6TwmiMi/ADisqptNZ4kQ0QCGApijqkMAnALgymNUItIF3r/wewO4FEAHEbnHbKrws1u57wfQq8HlZN821xKRdvAWe76q/sN0HoOuBXCziOyFd7pujIi8bjaSUWUAylT17F9yS+AtezcaB+A7VS1X1Z8A/APANYYzhZ3dyn0TgD4i0ltEYuA9KPKO4UzGiPebhOcD2KWqz5vOY5Kq/klVk1U1Dd5/Fx+pquP3zpqiqocA7BORvr5NYwHsNBjJpFIAI0Qk3vc7MxYuOLgcbTpAS6hqjYj8HsByeI94L1DVHYZjmXQtgMkAtovINt+2/1DV9w1mosjxEIB8347QtwDuM5zHCFXdKCJLAGyB9x1mW+GCT6ryE6pERA5kt2kZIiIKAsudiMiBWO5ERA7EciciciCWOxGRA7HciYgciOVORORALHciIgf6fwF9eHsBA+oUAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# iterate for robust location\n",
    "mu, w = [], np.ones_like(x)\n",
    "ilist = range(10)\n",
    "for _ in ilist:\n",
    "    m = sum(w*x) / sum(w)\n",
    "    mu.append(m) # save for plotting\n",
    "    w = Wc(x-m)\n",
    "    print (w[:3].T)\n",
    "\n",
    "plt.plot(ilist, mu,'ko-'); \n",
    "plt.plot(ilist, x[1:].mean()*np.ones_like(mu), 'b-', alpha=0.5); \n",
    "plt.ylim(-0.2,2.2)\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dispersion\n",
    "\n",
    "- MLE looks for the parameter for the optimial likelihood function\n",
    "\n",
    ">$\\displaystyle L(\\sigma) = p(D \\lvert \\sigma) = \\prod_{i=1}^n \\ell_i(\\sigma)$\n",
    "\n",
    "> with $\\displaystyle\\quad\\ell_i(\\sigma) = \\frac{1}{\\sigma}\\ f_0\\!\\left(\\frac{x_i}{\\sigma}\\right)$\n",
    "\n",
    ">$\\displaystyle \\hat{\\sigma} = \\arg \\max_{\\sigma} \\frac{1}{\\sigma^n}\\prod_i\\ f_0\\!\\left(\\frac{x_i}{\\sigma}\\right)$\n",
    "\n",
    "- Very general but for \"scale\" family of functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Influence and $\\rho$ functions\n",
    "\n",
    "- Influence function\n",
    "\n",
    ">$\\displaystyle \\psi = -f_0'\\big/f_0$\n",
    "\n",
    "- The $\\rho$-function\n",
    "\n",
    ">$\\displaystyle \\rho(t) = t\\, \\psi(t)$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimum\n",
    "\n",
    "- Derivative vanishes at optimal $\\hat{\\sigma}$\n",
    "\n",
    ">$\\displaystyle \\frac{1}{n} \\sum_{i=1}^n \\rho\\left(\\frac{x_i}{\\hat{\\sigma}}\\right) = 1$\n",
    "\n",
    "- Often $\\rho$ is normalized s.t. it tends to 1 as the argument goes to $\\infty$\n",
    "\n",
    "> Scaling is accommodated by new parameter $\\delta$\n",
    "\n",
    ">$\\displaystyle \\frac{1}{n} \\sum_{i=1}^n \\rho\\left(\\frac{x_i}{\\hat{\\sigma}}\\right) = \\delta$\n",
    "\n",
    "> For example, $\\delta=1/2\\ $ is often used "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What $\\rho$?\n",
    "\n",
    "- For Gaussian likelihood\n",
    "\n",
    "\n",
    ">$\\displaystyle f_0(t) = \\beta\\, e^{-t^2/2} \\quad$ \n",
    "and\n",
    "\n",
    ">$\\displaystyle f_0'(t) = -t \\ \\beta\\, e^{-t^2/2}$ \n",
    "\n",
    "> hence\n",
    "\n",
    ">$\\psi(t) = t\\quad$ and $\\quad\\rho(t)=t^2$ \n",
    "\n",
    "> and we have\n",
    "\n",
    ">$\\displaystyle \\frac{1}{n} \\sum_{i=1}^n \\left(\\frac{x_i}{\\hat{\\sigma}}\\right)^2 = 1$\n",
    "\n",
    "- So $\\hat{\\sigma}$ is the **RMS** (as before)\n",
    "\n",
    ">$\\displaystyle \\hat{\\sigma}^2 = \\frac{1}{n} \\sum_{i=1}^n x_i^2 $\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Weight\n",
    "\n",
    "- Define a new function $W$ \n",
    "\n",
    ">$ W(t) = \\left\\{ \\begin{array}{ll}\n",
    "         \\rho(t) \\big/ t^2 & \\mbox{if $t\\neq{}0$}\\\\\n",
    "         \\rho''(0) & \\mbox{if $t=0$}\\end{array} \\right.  $\n",
    "         \n",
    "> so we have\n",
    "\n",
    ">$\\displaystyle \\hat{\\sigma}^2 = \\frac{1}{n\\delta} \\sum_{i=1}^n W\\!\\left(\\frac{x_i}{\\hat{\\sigma}}\\right)\\,x_i^2$\n",
    "\n",
    "> where $\\delta$ is a new parameter\n",
    "\n",
    "- Play the same iterative trick"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xl8lOW99/HPL5M9ZCEbZCVhDwgJEHYXVNzKElvRgnWrVtQeq21t+2htezye06etdvE8Fc9zsGptUQGtS1QU9RQ3FEjYTTAYQiAbISGEJJA91/NHok9KgUySmdwz9/zerxevV2bmytzf0cw3d+657+sSYwxKKaXsxc/qAEoppVxPy10ppWxIy10ppWxIy10ppWxIy10ppWxIy10ppWzIqXIXkStFpEhEikXk/rOMuU5ECkWkQESed21MpZRS/SF9necuIg5gP3AZUA7kASuMMYW9xowD1gOXGGOOi0i8Meao+2IrpZQ6F2f23GcBxcaYEmNMG7AWyDltzO3AKmPMcQAtdqWUspa/E2OSgLJet8uB2aeNGQ8gIpsBB/CQMebt059IRFYCKwHCwsJmTJw4cSCZlVLKZ23fvr3WGBPX1zhnyt0Z/sA4YAGQDHwoIlOMMfW9BxljVgOrAbKzs01+fr6LNq+UUr5BRA45M86ZwzIVQEqv28k99/VWDuQaY9qNMQfpPkY/zpkASimlXM+Zcs8DxolIuogEAsuB3NPGvEr3XjsiEkv3YZoSF+ZUSinVD32WuzGmA7gb2AjsA9YbYwpE5GERWdozbCNwTEQKgU3Aj40xx9wVWiml1Ln1eSqku+gxd6WU6j8R2W6Mye5rnF6hqpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNqTlrpRSNuRUuYvIlSJSJCLFInL/GR6/RURqRGRXz7/vuD6qUkopZ/n3NUBEHMAq4DKgHMgTkVxjTOFpQ9cZY+52Q0allFL95Mye+yyg2BhTYoxpA9YCOe6NpZQ1mts6qaxvpqW90+ooSg1Kn3vuQBJQ1ut2OTD7DOOuEZELgf3AD4wxZWcYo5RHOnTsJI+8XcS7hdW0dXYR5O/HVeeN5MdXTiQpKsTqeEr1mzPl7ozXgReMMa0icgfwLHDJ6YNEZCWwEiA1NdVFm1ZqcN4vOsrdz++kyxi+NSeVCSPCKaxqYF1eGZuKavjvG2cwZ3SM1TGV6hcxxpx7gMhc4CFjzBU9tx8AMMb86izjHUCdMSbyXM+bnZ1t8vPzBxRaKVf55EAttzydx9j4Yfzp5mwSe+2lHzp2ktuezafieDPP3z6baanDLUyqVDcR2W6Mye5rnDPH3POAcSKSLiKBwHIg97SNJfS6uRTY15+wSlmhor6Zu9bsYFRMKM/fPvsfih1gVEwYL9w+h7jwIG7/y3ZqGlstSqpU//VZ7saYDuBuYCPdpb3eGFMgIg+LyNKeYfeISIGI7AbuAW5xV2ClXKGry/D9tTvp7DI8eVM2UaGBZxwXFx7E6ptm0NjSzo9e3E1ff+kq5SmcOuZujNkAbDjtvl/0+voB4AHXRlPKfV7IO0xe6XEeXTaVtNiwc46dODKCB66ayEOvF5K7u5KcrKQhSqnUwOkVqsrnnGhu55G3i5g3JoZlM5Kd+p4b56aRmRLFf7y5j+Y2PU1SeT4td+VznvywhBPN7Ty4KAMRcep7HH7CzxZlUNPYyrOflro1n1KuoOWufEpNYytPbz7I4qkJTE485wld/2RmWjQXT4jjv94/QENLu5sSKuUaWu7KpzzxfjGtHV388LLxA/r++y6fwInmdv70YYmLkynlWlruymfUNrXy3NbDfGNaEqPjhg3oOc5LimTRlASe3lxKo+69Kw+m5a58xnNbDtPW0cUdF40Z1PPccdFomlo7WJ9f7qJkSrmelrvyCa0dnfx1yyEWTIhjbPzA9tq/NDU5iplpw/nzJwfp7NLz3pVn0nJXPuH13VXUNrVy6/x0lzzfrfPTKatr5t3Capc8n1KupuWubM8Yw9MfH2T8iGFcMC7WJc95+eSRJA8P4emPD7rk+ZRyNS13ZXs7y+oprGrg5nlpTp/X3heHn3DT3FFsK61jf3WjS55TKVfScle2t25bGSEBDpZmJrr0ea+ZnkyAQ1i7TZcuUJ5Hy13ZWlNrB6/vqWRJZgLhwQEufe6YYUFcPmkkL+8sp7VDpyRQnkXLXdla7q5KTrV1snyWexaHWT4rhfpT7Wws0A9WlWfRcle2tjbvMBNGhDMtJcotzz9/TCzJw0NYu+2wW55fqYHScle2VXy0kT3lJ7huZorLPkg9nZ+fcF12Cp8cOEb58VNu2YZSA6Hlrmwrd1clfgJLMhP6HjwIX5/WPb/7a7sq3bodpfpDy13ZkjGG13ZXMm9MLPHhwW7dVkp0KDPThvPyjnJdqUl5DC13ZUt7yk9w6Ngpl5/+eDZXT0viQM1JPqtoGJLtKdUXLXdlS6/tqiTQ4ccV540cku0tnpJIoMOPV3ZWDMn2lOqLlruync4uwxt7KlkwIY7IENee2342kaEBXDIxntzdlXR0dg3JNpU6Fy13ZTtbS45xtLF1yBeyvnpaErVNrXxcXDuk21XqTLTcle28tquSsEAHl2bED+l2L57Y/ZeCHppRnkDLXdlKa0cnb31WxRWTRxIc4BjSbQf5O1g0NYGNBUc42doxpNtW6nRa7spWPiiqoaGlgyVZQ3OWzOmuzkqipb1L53lXltNyV7aSu7uS6LBAzh/rmnnb+yt71HASI4PJ3a0XNClrabkr2zjZ2sF7+6pZNCWBAIc1P9p+fsKSrEQ+3F9D3ck2SzIoBVruykbeLaympb2LpRYdkvnS0sxEOroMG/ZWWZpD+TYtd2Ubr+2qICkqhBmpwy3NMSkhgrHxw/TQjLKUU+UuIleKSJGIFIvI/ecYd42IGBHJdl1EpfpWd7KNj76oZXFmAn5+7pkB0lkiQk5mItsO1lFZ32xpFuW7+ix3EXEAq4CrgEnAChGZdIZx4cC9wFZXh1SqLxv2VtHRZcjJHNoLl85mSc+cNq/r3ruyiDN77rOAYmNMiTGmDVgL5Jxh3L8DvwFaXJhPKafk7q5kbPwwMhLCrY4CQFpsGJkpUXpoRlnGmXJPAnqvAFzec99XRGQ6kGKMefNcTyQiK0UkX0Tya2pq+h1WqTOprG9m28E6cjIT3bYox0DkZCZSUNlA8dFGq6MoHzToD1RFxA/4PXBfX2ONMauNMdnGmOy4uLjBblopAN7Y0713vGSIpvd11uKpCfhJ96IhSg01Z8q9AkjpdTu5574vhQPnAe+LSCkwB8jVD1XVUHltVyWZKVGkxYZZHeUfxEcEM3dMDLm7K3URDzXknCn3PGCciKSLSCCwHMj98kFjzAljTKwxJs0YkwZsAZYaY/LdklipXoqPNlFQ2TBki3L0V05mEqXHTrGn/ITVUZSP6bPcjTEdwN3ARmAfsN4YUyAiD4vIUncHVOpccndXIgJLprp3ndSBuuK8kQQ6/HR9VTXk/J0ZZIzZAGw47b5fnGXsgsHHUqpvxhhe313J3NExxEe4d53UgYoMCWDBhDje2FPJg4sycFh8Dr7yHXqFqvJaeytOcLD2JDkWTzfQl5ysJI42trK15JjVUZQP0XJXXit3VyUBDuHKyZ55SOZLl2bEExbo0EMzakhpuSuv1NlleH1PJQsmxBMZOjTrpA5UcICDKyaP5K3Pqmjt6LQ6jvIRWu7KK207WEd1Q6vHniVzuqVZiTS0dPBBkV68p4aGlrvySrm7KwgNdLAwY4TVUZwyf2ws0WGBvKbTEaghouWuvE5bRxcb9h7h8kkjCAkc2nVSByrA4ceiKQm8V1hNk66vqoaAlrvyOh/sr+FEc7vli3L0V05WIq0dXbxbeMTqKMoHaLkrr/PqzgpiwgK5YJx3zU80PXU4SVEhetaMGhJa7sqrnGhu59191SzJTLRsndSB8vMTlmQm8tEXtRxrarU6jrI573p3KJ/39mdVtHV0cfU0z1iUo79yshLp7DJs+EwPzSj30nJXXuWVnRWMjg0jMznS6igDMnFkOOPih5G7q6LvwUoNgpa78hoV9c1sKanj6mlJHrUoR3+ICDlZieSVHqdC11dVbqTlrrzGl4teXJ3lnYdkvrS0Z51XXV9VuZOWu/IKxhhe2VnOjFHDSY0JtTrOoKTGhJKVEqVnzSi30nJXXqGwqoH91U1e+0Hq6XKyEtlX1cAX1bq+qnIPLXflFV7dWUGAQ1g8xbNngHTWoi/XV9VDM8pNtNyVx+vsMry2q3sGyOFhgVbHcYn48GDmjYnltV26vqpyDy135fE2F9dytLHV6z9IPd3SrEQO151ix+F6q6MoG9JyVx5vfX4ZUaEBLJwUb3UUl7rqvJGEBDh4aXuZ1VGUDWm5K49Wf6qNdwqquToriSB/75gB0lnhwQEsmprA67urONWmM0Uq19JyVx7t1Z0VtHV2cV12itVR3OK67BSaWjvYsFenI1CupeWuPNr6/HLOS4pgUmKE1VHcYmbacNJjw1ifr4dmlGtpuSuP9VnFCQqrGvimTffaoXs6gmuzk9l2sI6DtSetjqNsRMtdeaz1+WUE+vt9dbm+XS2bnozDT3hR996VC2m5K4/U0t7JqzsruHLySCJDA6yO41bxEcEsGB/H33aU09HZZXUcZRNa7sojvf3ZERpaOmz7QerprpuZQnVDK+8X1VgdRdmElrvySGu2HCItJpR5Y2KsjjIkLpkYT3x4EGu2HrI6irIJp8pdRK4UkSIRKRaR+8/w+J0isldEdonIxyIyyfVRla/4/EgD+YeO863Zo/Dz88552/srwOHHilmpfLC/hrK6U1bHUTbQZ7mLiANYBVwFTAJWnKG8nzfGTDHGZAGPAL93eVLlM9ZsOUSgvx/LZiRbHWVIrZiVip8Iz209bHUUZQPO7LnPAoqNMSXGmDZgLZDTe4AxpqHXzTBAZ0JSA9LU2sErOypYPDXBNpOEOWtkZDALM+JZn19Ga0en1XGUl3Om3JOA3udolffc9w9E5F9E5ADde+73nOmJRGSliOSLSH5NjX5wpP7ZKzsrONnWyY1zRlkdxRI3zBlF3ck23tIrVtUguewDVWPMKmPMGOB/AT87y5jVxphsY0x2XFycqzatbMIYw3NbDjE5MYKslCir41hi/phY0mJCWbNFP1hVg+NMuVcAvc9HS+6572zWAlcPJpTyTdsPHefzI43cMGeU1y6APVh+fsK3Zo8i/9Bx9lU19P0NSp2FM+WeB4wTkXQRCQSWA7m9B4jIuF43FwFfuC6i8hVPbz5IZEgAOVmJVkex1LXZyYQEOHhm80Groygv1me5G2M6gLuBjcA+YL0xpkBEHhaRpT3D7haRAhHZBfwQuNltiZUtldWd4u3PjrBiViqhgf5Wx7FUVGgg18xI4tWdldQ0tlodR3kpp95FxpgNwIbT7vtFr6/vdXEu5WOe2VyKnwg3z/PND1JPd+v8dNZsOcxftxzih5eNtzqO8kJ6haqyXENLO+vzy1g8NYGEyBCr43iE0XHDWJgRz3NbDtHSrqdFqv7TcleWW59XRlNrB7edP9rqKB7l1vPTOXayjVd3nuv8BaXOTMtdWaqjs4tnNpcyKz2aKcmRVsfxKHNHxzApIYKnPj6IMXpdoOofLXdlqbcLjlBR38xt56dbHcXjiAi3nZ/OF0ebeH+/XvSn+kfLXVnGGMOqTQcYHRvGwowRVsfxSEsyE0mIDOa/Nh2wOoryMlruyjKbio6yr6qBuxaMweEjsz/2V6C/HysvHM220jq2HayzOo7yIlruyhLGGB7/ezFJUSFcPc3ey+gN1vKZqcSEBfL4pmKroygvouWuLPFpyTF2HK7nzotGE+DQH8NzCQl0cNsF6Xy4v4Y95fVWx1FeQt9VyhKrNhUTFx7EtT6yjN5g3ThnFBHB/qzSvXflJC13NeR2Hj7O5uJjrLxgNMEBDqvjeIXw4ABumZ/OxoJq9lc3Wh1HeQEtdzXkfv/ufoaHBnD97FSro3iVb89LY1iQP4+9t9/qKMoLaLmrIbWl5BgffVHLdxeMJSzItycI66/hYYHcen46G/Ye4bOKE1bHUR5Oy10NGWMMv91YxIiIIG6cqxOEDcR3LkgnKjSA375TZHUU5eG03NWQeb+ohvxDx/neJeP0WPsARQQHcNdFY3i/qEbPe1fnpOWuhkRXl+HRjUWkRodynZ4hMyg3zU0jLjyIRzd+rnPOqLPScldDYsNnVRRWNfD9heMI9Ncfu8EICXRwzyVjySs9rnPOqLPSd5lyu5b2Tn7z9udMGBFOTpZejeoK35yZyqiYUH61YR8dnV1Wx1EeSMtdud0zm0spq2vmZ4szdA4ZFwn09+OBqzLYX93EC3llVsdRHkjLXblVTWMrqzYVszAjngvGxVkdx1aumDyC2enR/OHd/Zxobrc6jvIwWu7KrX7/bhEt7Z389GsZVkexHRHh54sncfxUG4///Qur4ygPo+Wu3KawsoF1eWXcNDeN0XHDrI5jS+clRXLtjGT+/EkppbUnrY6jPIiWu3KLri7DQ7kFRIQEcO+l46yOY2s/umICgQ4/Hnq9QE+NVF/Rcldu8dL2craV1vHTqzKIDA2wOo6txYcH88PLJ/B+UQ1vfXbE6jjKQ2i5K5c71tTK/35rHzPThrNsRrLVcXzCzXNHMSkhgn97vYDGFv1wVWm5Kzf41Vuf09TSwS+/PgU/PfVxSPg7/Pjf35jC0cZWfveOzhqptNyVi20pOcZL28u5/cLRjB8RbnUcn5KVEsWNc0bxl09LdcUmpeWuXOdUWwf3/20PycNDuOcS/RDVCj+6YgKxw4L4yUt7aOvQK1d9mVPlLiJXikiRiBSLyP1nePyHIlIoIntE5H9EROdz9UGPvF1E6bFTPLJsKiGBOuujFSKCA/jVN6bw+ZFG/qjnvvu0PstdRBzAKuAqYBKwQkQmnTZsJ5BtjJkKvAQ84uqgyrN9euAYf/6klFvmpTFvTKzVcXzapRkjWDYjmSfeP8DuMj0846uc2XOfBRQbY0qMMW3AWiCn9wBjzCZjzKmem1sAPUXChzS1dvDjl3aTFhPKT66cYHUcBfx88STihgVx34u7aWnvtDqOsoAz5Z4E9J6ZqLznvrO5DXjrTA+IyEoRyReR/JoanarULn75ZiEV9c389tpMQgN16TxPEBkSwG+WTaX4aBO/01WbfJJLP1AVkRuAbODRMz1ujFltjMk2xmTHxekkUnbw5p4qXthWxh0XjiE7LdrqOKqXi8bHceOcUTz50UE2FR21Oo4aYs6UewXQe+mc5J77/oGILAQeBJYaY1pdE095srK6U9z/8h6yUqK47/LxVsdRZ/DgogwmjgznvvW7qW5osTqOGkLOlHseME5E0kUkEFgO5PYeICLTgP+mu9h1F8EHtHd28b0XdgLwxxXTCHDoWbWeKDjAwePXT6O5rZPvr91FZ5fOPeMr+nxHGmM6gLuBjcA+YL0xpkBEHhaRpT3DHgWGAS+KyC4RyT3L0ymbeHRjEbvK6vnNNVNJiQ61Oo46h7Hx4fxbzmQ+LTnG438vtjqOGiJOffpljNkAbDjtvl/0+nqhi3MpD/bGnkpWf1jCjXNG8bUpCVbHUU64dkYyWw4c47H/2c+U5AgumTjC6kjKzfRvadUv+6oa+PGLe8geNZyfLz79cgflqUSEX359ChkjI7h37S4O6tzvtqflrpxWf6qNlX/NJyLEnydumE6gv/74eJOQQAf/feMM/P2ElX/Jp6m1w+pIyo303amc0t7Zxd3P76T6RCv/dcMM4sODrY6kBiAlOpRV10+npPYkP1inH7DamZa76pMxhgdf2cvHxbX8x9fPY3rqcKsjqUGYNzaWny/K4N3Can755j6r4yg30csJVZ8e/3sx6/PLueeSsVyXndL3NyiPd8v8dA7XNfP05oOkRIfw7fnpVkdSLqblrs7plZ3l/O7d/Xx9WhI/uEwvVLKTBxdlUFF/ioffKCQxKoQrJo+0OpJyIT0so87q759X8+MX9zBndDS/uWYqIrqqkp04/ITHvjmNzOQovvfCTj45UGt1JOVCWu7qjD49cIy71uwgIyGC1Tdl65kxNhUS6OCZW2aSFhPK7c/ms/PwcasjKRfRd6z6J7vK6vnOs3mkRofy7K2ziAgOsDqScqPhYYGsuW02seFB3PJMHvuqGqyOpFxAy139gz3l9dz89DZihgWx5juziQ4LtDqSGgLxEcGsuW02IQEObvjTVj4/ogXv7bTc1Ve2HzrOt57cSkSIP899ZzYjIvRcdl+SEh3K87fPJsDhx4rVWyioPGF1JDUIWu4KgLzSOm56aisxwwJZt3KuTgbmo0bHDWPdHXMICXBw/ZNb2VuuBe+ttNwV/7Ovmhuf2sqIyGDW3TGXxKgQqyMpC42KCWPdHXMJD/ZnxZNb+KRYz6LxRlruPm7ttsPc/pd8xsWHs27lXD0Uo4DuQzQv3jmXxKhgbn5mG6/t+qf1eZSH03L3UcYYHntvP/e/vJcLxsWxduUc4sKDrI6lPEhCZAgv3jmPaanDuXftLp78sMTqSKoftNx9UFtHFw+8vJfH3vuCZTOS+dPN2YQF6cXK6p9FhgTwl1tnsWhKAr/csI9/e71AJxvzEvqO9jG1Ta18d80OtpXWcffFY7nv8vF65ak6p+AAB39cMY34iCCe2VxK8dEm/rhiGlGhepqsJ9M9dx+yt/wES/74MbvL6/nP5Vn86IoJWuzKKX5+wr8umcyvvzGFLSXHyFm1maIjjVbHUueg5e4jXtlZzrL/+wl+IvztrnnkZCVZHUl5oeWzUlm7ci6n2jr5+hObeWtvldWR1FloudvcqbYO7v/bHn6wbjfTUqPIvXs+5yVFWh1LebEZo4bzxvfOZ/yIcO56bgcP5RbQ0t5pdSx1Gj3mbmP7qhq4+/kdlNSe5LsLxvCDy8YT4NDf52rwRkQEs+6OOfz6rc95ZnMp2w7W8cfrpzEmbpjV0VQPfafbkDGGv3xaSs6qzTS2dLDmttn85MqJWuzKpYL8Hfzrksk8dXM2VSeaWfx/PmZ9XhnG6Nk0nkDf7TZTVneKG5/axi9eK2D+mBjeuvcC5o+NtTqWsrFLM0bw9vcvJCslip/8bQ/feTaf6oYWq2P5PLHqt2x2drbJz8+3ZNt21NVleG7bYX69oXtNzJ8uyuD6Wal6NowaMp1dhj9/UsqjGz8n0OHHzxdPYtmMZP0ZdDER2W6Mye5znJa79yupaeLBVz7j05JjXDAull99YwrJw3XiL2WN0tqT/OSlPWwrrWPBhDj+Pec8nYjOhbTcfcCptg5WbSpm9YclBAc4+NmiDK7LTtE9JWW5ri7DX7cc4jdvf05nl+FfLh7LygtHExzgsDqa19NytzFjDBsLqvn3NwqpqG/mmunJ3H/VRJ0bRnmcqhPN/Meb+3hzTxWjYkJ5aMlkLp4Yb3Usr+ZsuTv1gaqIXCkiRSJSLCL3n+HxC0Vkh4h0iMiygQRWztlVVs/y1Vu4c812woP9WX/HXH53XaYWu/JICZEhrLp+Omtum43DT/j2n/O46eltFFbqSk/u1ueeu4g4gP3AZUA5kAesMMYU9hqTBkQAPwJyjTEv9bVh3XPvn0PHTvLIxiLe3FNFTFgg9y4cx/WzUvHX0xuVl2jt6OQvnxzi8U3FNLS08/WsJH54+Xj9fKifnN1zd+YipllAsTGmpOeJ1wI5wFflbowp7Xmsa0Bp1VkdbWzhiU0HWLPlEAEOP+65dBwrLxzNMJ3FUXmZIH8Ht184muuyU3jig2Ke2VzKG3uruGVeGndeNEbX63UxZxoiCSjrdbscmO2eOOpLlfXNrP6whBe2Haajy/DNmSl8/9JxxOtiGsrLRYYG8MBVGdw8N43fv7ufJz8q4a+fHuKGOancfuFo4sP1Z9wVhnT3T0RWAisBUlNTh3LTXqOs7hRPvH+Al7aXYQx8Y3oS310wlrTYMKujKeVSiVEh/PbaTO68aDSrNh3gqY8P8uynh1gxM4U7Lhqjyz0OkjPlXgGk9Lqd3HNfvxljVgOrofuY+0Cew652ldXzzOaDvLGnCocI35yZwp0XjdHjkcr2xsaH84dvZnHvpeN44v1intt6mOe3HWZpZhK3np/G5ESd6G4gnCn3PGCciKTTXerLgevdmspHdHR2sbGgmqc3H2T7oeMMC/Lnlnlp3H7BaEZG6p+myrekxYbxyLJM7rl0HKs/LOHF/HL+tqOc2enR3Hp+OgszRuDw02s4nOXUee4i8jXgMcABPG2M+aWIPAzkG2NyRWQm8AowHGgBjhhjJp/rOX35bJnqhhZe2l7Oc1sOUXmihVExodwyL41lM5IJDw6wOp5SHuHEqXbW5R/m2U8OUVHfTEp0CDfMHsU1M5KJHea7p/7qRUwepqOzi01FNazLO8zfPz9Kl4F5Y2L49vx0LpkYr3skSp1FR2cX7xRW88zmg+SVHifAIVw2aQTLZ6Zy/thY/HzsvaPl7iEO1p7kpe1lvJhfztHGVuLCg7h2RjLXZafoh6RK9dMX1Y2szSvj5R3lHD/VTvLwEJbNSCYnK4l0H3k/ablbqLK+mTf3VJG7u5K9FSfwE7h4QjzLZ6WyYEKczquu1CC1dnSysaCatdsO82nJMYyBzORIcrKSWJyZYOvTKbXch9ixplY2fHaE13dVsq20DoCpyZEszUxk8dRE/YBUKTepOtHM67sreW1XJQWVDfgJzB8by+KpCSzMGEGMzY7Pa7kPgdLak7xTeIR3C6vZfug4XQbGxg9jaWYiSzITfebPRKU8RfHRRl7dWclruysoq2vGTyA7LZorJo/k8kkjbDH1sJa7G3R2GfaU1/PevmreKajmi6NNAGQkRHDZpBFcOXkkGQnhOuWuUhYzxlBQ2cA7BUfYWFBNUXUjAJMTI1iYMYIFE+KYmhzllScyaLm7SEV9Mx8jyrpnAAAHnUlEQVTtr+GjL2r5uLiWE83tOPyEWWnRXD55BAsz7LE3oJSdldaeZGPBETYWHGFnWT3GQFRoAOePjeWi8XFcND7Oa6b20HIfoIaWdvJL6/hwfy0ffVHDgZqTAIyMCObC8bFcMC6OC8bFEhWqkxwp5Y2On2zjo+JaPtxfwwf7a6hpbAVg4shwLhwfx5zR0WSnRRPhodecaLk76VhTK3mldWw9WMe2g3UUVjVgDAT5+zF7dAwXjuv+zT42fpgeblHKZowx7Ktq5MMvavigqIbth47T1tmFn8DkxEhmp0cze3QMs9KiiQz1jLLXcj+Dri5DSW0TOw/Xs7OsnryDdV8dNw/y92N66nBmj45mVlo000cN1yXBlPIxLe2d7Dh8nK0ldWwpOcbOsnraOroQgYkjI5ieGsW01OFMS40iPSbMkguotNyB2qZWdh2uZ1dZ97/d5fU0tnQAMCzInxmjhjMrPZo5o6OZkhRFoL+ef66U+v9a2jvZXVbPlpI68krr2F1WT2Nrd4dEhgSQlRLFtNQoslK6/w3F4VqfKndjDFUnWiisbKCwqoHCygY+qzxB+fFmABx+woQR4WT1/E+YlhLFmLhhPnfZslJqcLq6DAdqvvzr/zg7D9dTVN3IlzWaEh3CpIQIJidGMjkxgkmJEYyMCHbpIV3blnt7ZxcHapq6i/zLMq9qoP5UOwAikB4TRkZCBJkpkWSlDGdKUiQhgXqIRSnlek2tHewp7z46UNDTSwdrT371eHRYYHfRJ3SXfUZCBOmxYQO+Ut2Vy+x5lFWbinnsvS+A7uPkExMiuOq8BCb1/MebODKcMF2CTik1RIYF+TNvTCzzxsR+dV9TawefVzVQUNlAQeUJCqsaeGZzKW2d3SuR/mxRBt+5YLRbc3ldC35tSgLpsWFMTowgLSZMF4hWSnmcYUH+ZKd1n1L5pbaOLoqPNlFU3UBmcpTbM3hduY8fEc74EeFWx1BKqX4J9PfrPsKQGDEk29PdXqWUsiEtd6WUsiEtd6WUsiEtd6WUsiEtd6WUsiEtd6WUsiEtd6WUsiEtd6WUsiEtd6WUsiEtd6WUsiEtd6WUsiEtd6WUsiEtd6WUsiGnyl1ErhSRIhEpFpH7z/B4kIis63l8q4ikuTqoUkop5/VZ7iLiAFYBVwGTgBUiMum0YbcBx40xY4E/AL9xdVCllFLOc2bPfRZQbIwpMca0AWuBnNPG5ADP9nz9EnCpuHLRQKWUUv3izGIdSUBZr9vlwOyzjTHGdIjICSAGqO09SERWAit7bjaJSNFAQgOxpz+3F9PX4nns8jpAX4unGsxrGeXMoCFdickYsxpYPdjnEZF8ZxaI9Qb6WjyPXV4H6GvxVEPxWpw5LFMBpPS6ndxz3xnHiIg/EAkcc0VApZRS/edMuecB40QkXUQCgeVA7mljcoGbe75eBvzdGGNcF1MppVR/9HlYpucY+t3ARsABPG2MKRCRh4F8Y0wu8BTwVxEpBuro/gXgToM+tONB9LV4Hru8DtDX4qnc/lpEd7CVUsp+9ApVpZSyIS13pZSyIa8udxH5noh8LiIFIvKI1XkGS0TuExEjIrFWZxkIEXm05//HHhF5RUSirM7UX31NteEtRCRFRDaJSGHP++NeqzMNhog4RGSniLxhdZbBEJEoEXmp532yT0TmumtbXlvuInIx3VfGZhpjJgO/tTjSoIhICnA5cNjqLIPwLnCeMWYqsB94wOI8/eLkVBveogO4zxgzCZgD/IsXvxaAe4F9Vodwgf8E3jbGTAQyceNr8tpyB+4Cfm2MaQUwxhy1OM9g/QH4CeC1n3AbY94xxnT03NxC9zUR3sSZqTa8gjGmyhizo+frRrpLJMnaVAMjIsnAIuBPVmcZDBGJBC6k++xCjDFtxph6d23Pm8t9PHBBzyyUH4jITKsDDZSI5AAVxpjdVmdxoVuBt6wO0U9nmmrDKwuxt55ZWqcBW61NMmCP0b3j02V1kEFKB2qAZ3oOMf1JRMLctbEhnX6gv0TkPWDkGR56kO7s0XT/yTkTWC8ioz314qk+XstP6T4k4/HO9TqMMa/1jHmQ7sMCzw1lNvXPRGQY8Dfg+8aYBqvz9JeILAaOGmO2i8gCq/MMkj8wHfieMWariPwncD/wc3dtzGMZYxae7TERuQt4uafMt4lIF92T8dQMVb7+ONtrEZEpdP9G390zkWYysENEZhljjgxhRKec6/8JgIjcAiwGLvXUX7Tn4MxUG15DRALoLvbnjDEvW51ngOYDS0Xka0AwECEia4wxN1icayDKgXJjzJd/Qb1Ed7m7hTcflnkVuBhARMYDgXjhjHHGmL3GmHhjTJoxJo3uH4DpnljsfRGRK+n+83mpMeaU1XkGwJmpNrxCz5TbTwH7jDG/tzrPQBljHjDGJPe8N5bTPbWJNxY7Pe/pMhGZ0HPXpUChu7bn0XvufXgaeFpEPgPagJu9cE/Rbh4HgoB3e/4K2WKMudPaSM4721QbFscaqPnAjcBeEdnVc99PjTEbLMyk4HvAcz07DyXAt921IZ1+QCmlbMibD8sopZQ6Cy13pZSyIS13pZSyIS13pZSyIS13pZSyIS13pZSyIS13pZSyof8HuD1KoYhzfFsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def rho(t,c=1.0):\n",
    "    return c*c * 0.5 * np.log(1 + np.square(t/c))\n",
    "\n",
    "def W(t,c=1.0):\n",
    "    return rho(t,c) / np.square(t)\n",
    "\n",
    "t = np.linspace(-6,6,300); \n",
    "plt.plot(t,W(t)); plt.ylim(0,0.6);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20.02061035643945, 0.9127818753073593)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "naive, cheat = np.sqrt(np.square(x).mean()), np.sqrt(np.square(x[1:]).mean())\n",
    "naive, cheat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20.02061035643945, 0.9578626069648577)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAHwVJREFUeJzt3Xt0VPXd7/H3NwS5qhCJ3JQEKVCjgpE9XmqP2ovIRdTWVmWB9TlHV6zVc/QRddVDl7raRZfHqu1jn64qVaqtqY89Xo4iFLU+trbrWDWkRFDwyEUgCBJEgTZoTPI9f2QHhzCBycwkOzP781pr1uz9m9/s/R1lfWbnN/u3t7k7IiISH0VRFyAiIj1LwS8iEjMKfhGRmFHwi4jEjIJfRCRmFPwiIjGj4BcRiRkFv4hIzCj4RURipjjqAlIZNmyYl5eXR12GiEjeWL58+Q53L02nb68M/vLycmpqaqIuQ0Qkb5jZxnT7aqhHRCRmFPwiIjGj4BcRiRkFv4hIzCj4RURi5pDBb2bHmtnLZva2mb1lZteH7SVm9qKZvRs+D+3k/VeEfd41syty/QHaVVdXU15eTlFREeXl5VRXV3fXrkRE8lo6R/zNwDx3rwBOB641swrg+8BL7j4eeClc34+ZlQC3A6cBpwK3d/YFkY3q6mqqqqrYuHEj7s7GjRupqqpS+IuIpHDI4Hf3re5eGy7vAVYDo4ELgUfCbo8AF6V4+3nAi+6+090/Al4EpuWi8GTz58+nsbFxv7bGxkbmz5+f612JiOS9Lo3xm1k5UAm8Bgx3963hS9uA4SneMhrYnLReH7al2naVmdWYWU1DQ0NXymLTpk1dahcRibO0g9/MBgNPAje4++7k17ztju1Z3bXd3Re6e+DuQWlpWrOO9xkzZkyX2kVE4iyt4DezvrSFfrW7PxU2f2BmI8PXRwLbU7x1C3Bs0voxYVtOLViwgIEDB+7XNnDgQBYsWJDrXYmI5L10zuox4CFgtbvfm/TSs0D7WTpXAM+kePvzwFQzGxr+qDs1bMupOXPmsHDhQtr/Uhg+fDgLFy5kzpw5ud6ViEjeS+eI/0zgcuCrZrYifMwA7gTONbN3ga+H65hZYGYPArj7TuBHwBvh44dhW87NmTOHv//97wDceuutCn0RkU4c8uqc7v5XwDp5+Wsp+tcAVyWtLwIWZVpgV4wePZqRI0fqyp4iIgdRcDN3gyDgjTfeiLoMEZFeq+CCP5FI8M4777B79+5DdxYRiaGCC/4gCABYvnx5xJWIiPROBRv8GucXEUmt4IK/tLSUsrIyjfOLiHSi4IIf2sb5dcQvIpJaQQZ/EARs2LCBHTt2RF2KiEivU5DBn0gkAP3AKyKSSkEG/ymnnAKgcX4RkRQKMviHDBnChAkTNM4vIpJCQQY/aAaviEhnCjb4E4kE77//Pu+//37UpYiI9CoFG/yayCUiklrBBn9lZSVFRUUKfhGRDgo2+AcNGkRFRYXG+UVEOijY4IfPZ/C23RJYRESgwIM/CAJ27NjBxo0boy5FRKTXSOeeu4vMbLuZrUpqezzpNozvmdmKTt77npmtDPv1+GB7+wxejfOLiHwunSP+h4FpyQ3ufqm7n+zuJwNPAk8d5P1fCfsGmZeZmUmTJtG3b1+N84uIJEnnnruvmFl5qtfMzIBLgK/mtqzc6NevH5MmTdIRv4hIkmzH+P8L8IG7v9vJ6w68YGbLzazqYBsysyozqzGzmoaGhizL+lwQBNTU1NDa2pqzbYqI5LNsg3828NhBXv+yu58CTAeuNbOzOuvo7gvdPXD3oLS0NMuyPpdIJNi9ezdr167N2TZFRPJZxsFvZsXAN4HHO+vj7lvC5+3A08Cpme4vU+0zeDXOLyLSJpsj/q8Da9y9PtWLZjbIzA5vXwamAqtS9e1OJ5xwAv3799c4v4hIKJ3TOR8DXgUmmlm9mV0ZvnQZHYZ5zGyUmS0NV4cDfzWzOuB1YIm7L8td6ekpLi6msrJSwS8iEkrnrJ7ZnbT/S4q294EZ4fJ6YHKW9eVEIpHgwQcfpLm5meLiQ35kEZGCVtAzd9sFQUBjYyNr1qyJuhQRkcjFIvjbZ/DqB14RkZgE/4QJEzj88MM1zi8iQkyCv6ioiClTpuiIX0SEmAQ/tI3z19XV0dTUFHUpIiKRik3wJxIJmpqaWLlyZdSliIhEKjbBr3vwioi0iU3wjx07lpKSEo3zi0jsxSb4zWzflTpFROIsNsEPbeP8q1atorGxMepSREQiE6vgD4KAlpYW6urqoi5FRCQysQp+zeAVEYlZ8I8aNYoRI0ZonF9EYi1WwW9mJBIJHfGLSKzFKvihbZz/nXfeYffu3VGXIiISidgFfyKRwN2pra2NuhQRkUjELvg1g1dE4i6dWy8uMrPtZrYqqe0OM9tiZivCx4xO3jvNzN4xs7Vm9v1cFp6p0tJSysrKNM4vIrGVzhH/w8C0FO0/dfeTw8fSji+aWR/gF8B0oAKYbWYV2RSbK5rBKyJxdsjgd/dXgJ0ZbPtUYK27r3f3JuA/gAsz2E7OJRIJ1q9fz4cffhh1KSIiPS6bMf7rzOzNcChoaIrXRwObk9brw7aUzKzKzGrMrKahoSGLsg6tfZx/+fLl3bofEZHeKNPg/yUwDjgZ2Arck20h7r7Q3QN3D0pLS7Pd3EFNmTIF0AxeEYmnjILf3T9w9xZ3bwV+RduwTkdbgGOT1o8J2yI3ZMgQxo8fr3F+EYmljILfzEYmrX4DWJWi2xvAeDMba2aHAZcBz2ayv+6gGbwiElfpnM75GPAqMNHM6s3sSuAuM1tpZm8CXwH+New7ysyWArh7M3Ad8DywGvi9u7/VTZ+jy4IgYMuWLWzdujXqUkREelTxoTq4++wUzQ910vd9YEbS+lLggFM9e4P2K3XW1NQwa9asiKsREek5sZu5266yspKioiKN84tI7MQ2+AcNGkRFRYXG+UUkdmIb/PD5DF53j7oUEZEeE+vgTyQSNDQ0sGnTpqhLERHpMbEOfl2pU0TiKNbBP3nyZPr27atxfhGJlVgHf79+/TjppJN0xC8isRLr4Ie2cf6amhpaW1ujLkVEpEfEPviDIGDXrl2sW7cu6lJERHpE7IO/fQavxvlFJC5iH/wVFRX0799f4/wiEhuxD/6+fftSWVmpI34RiY3YBz+0jfPX1tbS0tISdSkiIt1OwU/bOH9jYyOrV6+OuhQRkW6n4EczeEUkXhT8wMSJExk8eLDG+UUkFtK5A9ciM9tuZquS2n5iZmvM7E0ze9rMhnTy3vfCO3WtMLNeezhdVFTElClTdMQvIrGQzhH/w8C0Dm0vAie6+yTg/wG3HuT9X3H3k909yKzEnpFIJFixYgVNTU1RlyIi0q0OGfzu/gqws0PbC+E9dQH+BhzTDbX1qCAIaGpqYtWqVPeNFxEpHLkY4/9vwB86ec2BF8xsuZlV5WBf3UYzeEUkLrIKfjObDzQD1Z10+bK7nwJMB641s7MOsq0qM6sxs5qGhoZsysrI2LFjKSkp0Ti/iBS8jIPfzP4FOB+Y453cu9Ddt4TP24GngVM72567L3T3wN2D0tLSTMvKmJkRBIGO+EWk4GUU/GY2DbgFuMDdGzvpM8jMDm9fBqYCvXoAPZFIsGrVKvbu3Rt1KSIi3Sad0zkfA14FJppZvZldCfw7cDjwYniq5v1h31FmtjR863Dgr2ZWB7wOLHH3Zd3yKXIkCAJaWlpYsWJF1KWIiHSb4kN1cPfZKZof6qTv+8CMcHk9MDmr6npY+w+8NTU1nHHGGRFXIyLSPTRzN8moUaMYMWKExvlFpKAp+JOY2b5bMYqIFCoFfwdBELBmzRr27NkTdSkiIt1Cwd9BIpHA3amtrY26FBGRbqHg76D9Es0a5xeRQqXg76C0tJSysjKN84tIwVLwp6AZvCJSyBT8KSQSCdavX8/OnTsP3VlEJM8o+FPQrRhFpJAp+FOYMmUKoOAXkcKk4E9hyJAhjB8/XuP8IlKQFPyd0AxeESlUCv5OBEFAfX0927Zti7oUEZGcUvB3IvlKnSIihUTB34nKykqKioo0zi8iBUfB34lBgwZRUVGhI34RKTgK/oNon8HbyS2FRUTyUlrBb2aLzGy7ma1KaisxsxfN7N3weWgn770i7POumV2Rq8J7QiKRoKGhgc2bN0ddiohIzqR7xP8wMK1D2/eBl9x9PPBSuL4fMysBbgdOA04Fbu/sC6I30pU6RaQQpRX87v4K0PHCNRcCj4TLjwAXpXjrecCL7r7T3T8CXuTAL5Bea/LkyfTt21fj/CJSULIZ4x/u7lvD5W3A8BR9RgPJ4yT1YVte6NevHyeddJKO+EWkoOTkx11v+/Uzq19AzazKzGrMrKahoSEXZeVE+wxe/cArIoUim+D/wMxGAoTP21P02QIcm7R+TNh2AHdf6O6BuwelpaVZlJVbQRCwa9cu1q5dG3UpIiI5kU3wPwu0n6VzBfBMij7PA1PNbGj4o+7UsC1vaAaviBSadE/nfAx4FZhoZvVmdiVwJ3Cumb0LfD1cx8wCM3sQwN13Aj8C3ggfPwzb8kZFRQX9+/fXOL+IFAzrjWPXQRB4bzrC/tKXvkRxcTGvvPJK1KWIiKRkZsvdPUinr2bupiEIAmpra2lpaYm6FBGRrCn405BIJPjnP//JmjVroi5FRCRrCv40aAaviBQSBX8aJk6cyODBg3Vmj4gUBAV/GoqKipgyZYqO+EWkICj405RIJKirq6OpqSnqUkREsqLgT1MQBHz66aesWrXq0J1FRHoxBX+aNINXRAqFgj9NY8eOpaSkROP8IpL3FPxpMjOCINARv4jkPQV/FwRBwMqVK9m7d2/UpYiIZEzB3wWJRIKWlhbq6uqiLkVEJGMK/i7QDF4RKQQK/i4YPXo0I0aM0Di/iOQ1BX8XtP/AqyN+EclnCv4uSiQSrFmzhj179kRdiohIRhT8XRQEAe5ObW1t1KWIiGQk4+A3s4lmtiLpsdvMbujQ5xwz25XU57bsS45W+w+8GucXkXxVnOkb3f0d4GQAM+sDbAGeTtH1L+5+fqb76W2OPvpoxowZo3F+EclbuRrq+Rqwzt035mh7vVoikdARv4jkrVwF/2XAY528doaZ1ZnZH8zshBztL1JBELBu3Tp27twZdSkiIl2WdfCb2WHABcD/TvFyLVDm7pOBnwP/5yDbqTKzGjOraWhoyLasbtV+pc7ly5dHXImISNfl4oh/OlDr7h90fMHdd7v7P8LlpUBfMxuWaiPuvtDdA3cPSktLc1BW95kyZQqgGbwikp9yEfyz6WSYx8xGmJmFy6eG+/swB/uM1JAhQxg/frzG+UUkL2V8Vg+AmQ0CzgWuTmr7LoC73w98C7jGzJqBvcBl7u7Z7LO3CIKAv/zlL1GXISLSZVkd8bv7P939KHffldR2fxj6uPu/u/sJ7j7Z3U939/+bbcG9RSKRoL6+nm3btkVdiohIl2jmboY0kUtE8pWCP0OVlZUUFRUp+EUk7yj4MzR48GCOP/54ndkjInlHwZ+F9hm8BfJ7tYjEhII/C0EQsH37djZv3hx1KSIiaVPwZ6F9Bq/G+UUknyj4szBp0iSKi4s1zi8ieUXBn4X+/fszadIkHfGLSF5R8GcpCAL9wCsieUXBn6VEIsHHH3/MunXroi5FRCQtCv4saQaviOQbBX+W3nzzTQBmz55NeXk51dXVEVckInJwCv4sVFdXc8011+xb37hxI1VVVQp/EenVFPxZmD9/Po2Njfu1NTY2Mn/+/IgqEhE5NAV/FjZt2tSldhGR3kDBn4UxY8akbD/mmGN6uBIRkfQp+LOwYMECBg4ceED7yJEjaW1tjaAiEZFDyzr4zew9M1tpZivM7IBzGq3NfWa21szeNLNTst1nbzFnzhwWLlxIWVkZZkZZWRmzZ8/m9ddf59Zbb426PBGRlLK6526Sr7j7jk5emw6MDx+nAb8MnwvCnDlzmDNnzr51d+fII4/krrvuYty4cVRVVUVYnYjIgXIV/AdzIfCb8CbrfzOzIWY20t239sC+e5yZ8fOf/5yNGzfyve99jzFjxjBt2rSoyxIR2ScXY/wOvGBmy80s1eHtaCD5gvX1YVvBKi4u5vHHH+fEE0/kkksuoa6uLuqSRET2yUXwf9ndT6FtSOdaMzsrk42YWZWZ1ZhZTUNDQw7Kitbhhx/OkiVLOOKII5g5cyZbtmyJuiQRESAHwe/uW8Ln7cDTwKkdumwBjk1aPyZs67idhe4euHtQWlqabVm9wujRo1myZAm7du3i/PPPZ8+ePVGXJCKSXfCb2SAzO7x9GZgKrOrQ7VngO+HZPacDuwp1fD+VyZMn8/vf/56VK1dy2WWX0dzcHHVJIhJz2R7xDwf+amZ1wOvAEndfZmbfNbPvhn2WAuuBtcCvgO9luc+8M336dH7xi1+wdOlSrr/+el27X0QildVZPe6+Hpicov3+pGUHrs1mP4Xg6quvZt26dfzkJz9h3Lhx3HjjjVGXJCIx1ROnc0rozjvvZMOGDdx0002UlZVx8cUXR12SiMSQLtnQg4qKivjNb37Daaedxty5c3nttdeiLklEYkjB38MGDBjAM888w6hRo5g1axYbNmyIuiQRiRkFfwSOPvpoli5dSnNzMzNmzOCjjz6KuiQRiREFf0QmTpzI008/zbp16/jmN79JU1NT1CWJSEwo+CN09tlns2jRIv70pz9x1VVX6TRPEekROqsnYnPnzmXDhg3cdtttjBs3jttvvz3qkkSkwCn4e4Ef/OAHrF+/njvuuIPjjjuOyy+/POqSRKSAKfh7ATPjgQceYNOmTVx55ZUce+yxnHPOOVGXJSIFSmP8vcRhhx3Gk08+yRe+8AW+8Y1vsHr16qhLEpECpeDvRYYMGcLSpUs57LDDmDlzJtu3b4+6JBEpQAr+Xqa8vJznnnuObdu2ccEFF7B3796oSxKRAqPg74USiQS/+93veP3115k7dy6tra1RlyQiBUTB30tddNFF3HPPPTz11FPccsstUZcjIgVEZ/X0YjfccAPr16/nnnvuYdy4cVxzzTVRlyQiBUDB34uZGT/72c947733uO666ygrK2PGjBlRlyUieU5DPb1cnz59eOyxx5g8eTKXXnopK1asiLokEclzGQe/mR1rZi+b2dtm9paZXZ+izzlmtsvMVoSP27IrN54GDx7Mc889x9ChQ5k5cyb19fVRlyQieSybI/5mYJ67VwCnA9eaWUWKfn9x95PDxw+z2F+sjRo1iiVLlrBnzx5mzpzJ7t27oy5JRPJUxsHv7lvdvTZc3gOsBkbnqjA50EknncQTTzzBW2+9xaWXXkpzc3PUJYlIHsrJGL+ZlQOVQKp7CZ5hZnVm9gczO+Eg26gysxozq2loaMhFWQVp6tSp3H///SxbtozzzjuPsrIyioqKKC8vp7q6OuryRCQPWLbXgDezwcCfgQXu/lSH144AWt39H2Y2A/g3dx9/qG0GQeA1NTVZ1VXoLrjgAhYvXrxf28CBA1m4cCFz5syJqCoRiYqZLXf3IJ2+WR3xm1lf4EmgumPoA7j7bnf/R7i8FOhrZsOy2ae0qaurO6CtsbGR+fPnR1CNiOSTbM7qMeAhYLW739tJnxFhP8zs1HB/H2a6T/nc5s2bU7Zv3LiRe++9l3fffbeHKxKRfJHNEf+ZwOXAV5NO15xhZt81s++Gfb4FrDKzOuA+4DLX/QVzYsyYMSnb+/bty7x585gwYQJf/OIXufnmm3nllVf0Q7CI7JP1GH930Bj/oVVXV1NVVUVjY+O+tvYx/jPPPJPnnnuOxYsX8/LLL/PZZ59RUlLC9OnTmTVrFtOmTePII4+MsHoRybWujPEr+PNYdXU18+fPZ9OmTYwZM4YFCxYc8MPunj17eOGFF1i8eDFLlixhx44dFBcXc9ZZZzFr1ixmzZrFuHHjIvoEIpIrCn5JqaWlhb/97W8sXryYxYsX8/bbbwNw/PHH7/sSOOOMM+jTp0/ElYpIVyn4JS3r16/f9yXw5z//mebmZo466ihmzJjBrFmzOO+88zjiiCOiLlNE0tBjp3NKfjvuuOO4/vrr+eMf/8iOHTt4/PHHmT59OkuWLOGSSy5h2LBhnHvuudx3331s2LBh3/uqq6spLy/XxDGRPKUjfjlAc3Mzr7766r6/BtasWQPAiSeeyHHHHcfzzz/Pp59+uq+/Jo6JRE9DPZJTa9eu3fcl8PLLL6fsU1JSwq9+9SuGDx++7zF48GDCaRwi0s0U/NJtioqKSPffzIABAzj66KP3+zLouN7eVlJSktaXRDpnMonEUVeCH3fvdY8pU6a49E5lZWUOHPAYPXq019bW+rJly/yRRx7xu+66y+fNm+dz5871qVOn+uTJk33EiBHep0+flO8vLi72UaNGeWVlpU+bNs2vuOIKv/nmm/3uu+/23/72t/7CCy/4j3/8Yx8wYMB+7xs4cKA/+uij3f65H330US8rK3Mz87Kysh7ZZ5T7jXLf+syZ7Reo8TQzVkf80iUHmziWzpF3a2srH374Idu3b+eDDz7Y9+i43v5oamo65Db79OnDhAkTGDBgAP3796d///77LXdc7+ry0qVLmTdvHnv37s3oM2cq2//W+bhvfebM96uhHulWPTXc4u7s3r1735fA2Wef3ekw07e//W327t3LJ598wieffHLQ5Vz9mzczBg0aRHFxMX369Nn3nLycqu1Qr7c/L168eL8waDd48GC+853vUFRUhJml/dyVvgsWLOCjjz46YN9Dhw7lRz/60b7ttf93SF5O1Zbu8g033MCHHx54Oa9hw4Zx33337Tcc2HFoMJt1M+Pqq69mx44dB+y7tLSUBx544KDb6qgrr1911VWkuhR9WVkZ77333kG302GbMQ3+Zctg27bcFyS9wk033cSOVKFw1FHcfffdaW3D3WlpaeGzzz6jqamJzz77bL9Hclv78qJf/7rT7Z03dSotra14ayutra20tLTQ2tpKqzut7ctJj5aWlk5fa21t3betlpYWtn3wQaf7HTxoEND2FxRAqzsc6k/8tP4LSZS2Ac+Hy2a27/9vOroS/MUZ1CYSiYsvvpiHH36YT5OGf/oddhgXX3xx2tswM4qLiykuLmbAgAFpvefZZ5/t9Atn9uzZae+7q3LxRddR+5cAHPzLYv4PfsDOnTsPeH/J0KHccccd+/3V1L6csq1tJeXrqfrfeeedfPTxxwfsd+iQIdxyyy0H9O9sPdXn3m99/xcBuOeee/h4164D3jvkyCO58cYbD7r9bGr56U9/yq4Ut1Lt7EKMOZHujwE9+dCPu9KZKH58e/TRR33gwIE9/qNyVPuNct/6zJnvly78uBt5yKd6KPilt8nnsz3ybd/6zDqrR0RE0qBr9YiISKcU/CIiMZPtzdanmdk7ZrbWzL6f4vV+ZvZ4+PprZlaezf5ERCR72dxsvQ/wC2A6UAHMNrOKDt2uBD5y9y8APwX+V6b7ExGR3MjmiP9UYK27r3f3JuA/gAs79LkQeCRcfgL4mulyjSIikcom+EcDm5PW68O2lH3cvRnYBRyVamNmVmVmNWZWk2r6soiI5Eavmbnr7guBhQBm1mBmGzPc1DDgwAtuFDZ95sIXt88L+sxdVZZux2yCfwtwbNL6MWFbqj71ZlYMHAkcOAe9A3cvzbQoM6tJ91zWQqHPXPji9nlBn7k7ZTPU8wYw3szGmtlhwGXAsx36PAtcES5/C/hP740zxkREYiTjI353bzaz62i7mFwfYJG7v2VmP6Rt6vCzwEPAb81sLbCTti8HERGJUFZj/O6+FFjaoe22pOVPgG9ns48MLOzh/fUG+syFL26fF/SZu02vvFaPiIh0H12yQUQkZgom+A91+YhCY2bHmtnLZva2mb1lZtdHXVNPMbM+ZvZ3M3su6lp6gpkNMbMnzGyNma02szOirqm7mdm/hv+uV5nZY2bWP+qacs3MFpnZdjNbldRWYmYvmtm74fPQ7th3QQR/mpePKDTNwDx3rwBOB66NwWdudz2wOuoietC/Acvc/YvAZAr8s5vZaOB/AIG7n0jbySOFeGLIw8C0Dm3fB15y9/HAS+F6zhVE8JPe5SMKirtvdffacHkPbWHQceZ0wTGzY4CZwINR19ITzOxI4CzazpDD3Zvc/cB7ExaeYmBAOP9nIPB+xPXknLu/QtvZjsmSL3PzCHBRd+y7UII/nctHFKzwqqeVwGvRVtIjfgbcAqR/F+r8NhZoAH4dDm89aGaDoi6qO7n7FuBuYBOwFdjl7i9EW1WPGe7uW8PlbcDw7thJoQR/bJnZYOBJ4AZ3P/COzQXEzM4Htrv78qhr6UHFwCnAL929Evgn3fTnf28RjmtfSNuX3ihgkJnNjbaqnhdOdu2W0y4LJfjTuXxEwTGzvrSFfrW7PxV1PT3gTOACM3uPtuG8r5rZo9GW1O3qgXp3b/9r7gnavggK2deBDe7e4O6fAU8BX4q4pp7ygZmNBAift3fHTgol+NO5fERBCS9v/RCw2t3vjbqenuDut7r7Me5eTtv/4/9094I+EnT3bcBmM5sYNn0NeDvCknrCJuB0MxsY/jv/GgX+g3aS5MvcXAE80x076TVX58xGZ5ePiLis7nYmcDmw0sxWhG3/M5xNLYXlvwPV4UHNeuC/RlxPt3L318zsCaCWtrPX/k4BzuI1s8eAc4BhZlYP3A7cCfzezK4ENgKXdMu+NXNXRCReCmWoR0RE0qTgFxGJGQW/iEjMKPhFRGJGwS8iEjMKfhGRmFHwi4jEjIJfRCRm/j/sIQp6R3WVYwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "scale, delta = [np.sqrt(np.square(x).mean())], 0.3\n",
    "for _ in range(10):\n",
    "    w = W(x/scale[-1])\n",
    "    scale.append(np.sqrt((np.square(x)*w).mean()/delta))\n",
    "plt.plot(scale,'ko-'); \n",
    "plt.plot(cheat*np.ones_like(scale), 'r-', alpha=0.5); \n",
    "scale[0], scale[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### New opportunities\n",
    "\n",
    "\n",
    "- Robust regression\n",
    "\n",
    "> Minimize robust $M$-scale of residuals instead of RMS\n",
    "\n",
    "- Robust PCA\n",
    "\n",
    "> Maximize robust $M$-scale (squared) instead of variance (Maronna 2005)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
